setwd("D:\\R")

rt=read.table("GSE42546_CleanedRawCounts.txt",sep="\t",header=T,check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
rt=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
rt<-rt[!apply(rt,1,function(x){sum(floor(x)==0)>3}),]#有三个表达为0，即去除。
##提取分组信息
col_names <- colnames(rt)  
sorted_col_names<- sort(col_names) 
rt<-rt[, sorted_col_names]
rt<-rt[, 17:62] #选择CD两组分析，即"控制"与“抑郁”
group <- substring(colnames(rt),1,1)
group <- factor(group,labels = c("控制","抑郁"))


library(DESeq2)
coldata <- data.frame(row.names = colnames(rt), group)
dds <- DESeqDataSetFromMatrix(countData = rt, colData = coldata, design= ~group)
dds1<-DESeq(dds)
res <- results(dds1)
res <- res[order(res$log2FoldChange, res$padj, decreasing = c(FALSE, TRUE)), ]
res
##筛选差异基因
rownames(result)
res[which(res$log2FoldChange >= 1 & res$padj < 0.05),'sig'] <- 'up'
res[which(res$log2FoldChange <= -1 & res$padj < 0.05),'sig'] <- 'down'
res[which(abs(res$log2FoldChange) <= 1 | res$padj >= 0.05),'sig'] <- 'none'
res
res_select <- subset(res, sig %in% c('up', 'down'))
result<-data.frame(res_select)
result$GeneNames <- row.names(result)  
library(ggplot2)

#横轴展示log2FoldChange，纵轴展示 -log10 转化后的 padj
p <- ggplot(data = res, aes(x = log2FoldChange, y = -log10(padj), color = sig)) +
  geom_point(size = 1) +  #绘制散点图
  scale_color_manual(values = c('red', 'gray', 'green'), limits = c('up', 'none', 'down')) +  #自定义点的颜色
  labs(x = 'log2 Fold Change', y = '-log10 adjust p-value', title = '火山图', color = '') +  #坐标轴标题
  theme(plot.title = element_text(hjust = 0.5, size = 14), panel.grid = element_blank(), #背景色、网格线、图例等主题修改
        panel.background = element_rect(color = 'black', fill = 'transparent'), 
        legend.key = element_rect(fill = 'transparent')) +
  geom_vline(xintercept = c(-1, 1), lty = 3, color = 'black') +  #添加阈值线（x=1,y=1）
  geom_hline(yintercept = 1, lty = 3, color = 'black') +
  xlim(-5, 5) + ylim(0, 10)  
p
#GO富集分析
library(org.Hs.eg.db)
library(clusterProfiler)
library(ggplot2)
library(forcats)
ego <- enrichGO(gene = result$GeneNames,
                    OrgDb=org.Hs.eg.db,
                    keyType = "SYMBOL",
                    ont = "ALL",
                    pAdjustMethod = "BH",
                    minGSSize = 1,
                    pvalueCutoff = 0.05,
                    qvalueCutoff = 0.05,
                    readable = TRUE)
all <- as.data.frame(ego)
all$Description <- as.factor(all$Description) 
ggplot(all, aes(x = -log10(pvalue), y = fct_reorder(Description, -log10(pvalue), .fun = median, .na_rm = TRUE), group = ONTOLOGY)) +
  geom_point(aes(size = Count, fill = p.adjust), shape = 21, color = 'black', na.rm = TRUE) +
  facet_grid(ONTOLOGY ~ ., scale = 'free_y', space = 'free_y') +
  scale_fill_gradient(low = 'blue', high = 'red') +
  labs(title = 'GO Enrichment', y = 'GO term', x = '-log10(pvalue)') +
  guides(fill = guide_colorbar(reverse = TRUE, order = 1)) +
  theme_bw() +
  coord_cartesian(xlim = c(0, max(-log10(all$pvalue))))
ggplot(all, 
       aes(x=Description,y=Count, fill=ONTOLOGY)) + 
  geom_bar(stat="identity", width=0.8) +  
  scale_fill_manual(values = c("#6666FF", "#33CC33", "#FF6666") ) +  
  facet_grid(ONTOLOGY~., scale = 'free_y', space = 'free_y')+
  coord_flip() +  
  xlab("GO term") +  
  ylab("count") +  
  labs(title = "GO Terms Enrich")+  
  theme_bw()


